## "China's Ideological Spectrum" by Pan and Xu (2018)
## Replication by Chen & Chen 
# install.packages("psych")

## set working directory, where the ReadMe.txt is in## 


########### Show Results ###################

rm(list=ls(all=TRUE)) 
load("output/result_cfa_search_dropQ6.RData")
post.cov <- which(results.final$post.check == 1)
good.results <- results.final[post.cov,]
(best.id<-which.min(good.results$chisq))
which.max(good.results$cfi)
which.max(good.results$tli)
which.min(good.results$rmsea)
print(good.results[best.id, ])

## Assignment: 1, 2, 2, 3, 1, 2, 2
## ------------------------------------
f1 <- c(1, 2, 3, 4, 5, 10, 12, 13, 14, 17) # political liberalism
f5 <- c(24, 7, 8, 41, 44, 45, 50) # individual freedom
## ------------------------------------
f2 <- c(21, 25, 27, 29, 30, 37, 40) # free market
f3 <- c(32, 42, 43, 46, 47, 48, 49) # traditional values
f6 <- c(23, 26, 28, 35, 39) # national interest
f7 <- c(22, 31, 33, 34, 36, 38) # social justice
## ------------------------------------
f4 <- c(9, 11, 15, 16, 18, 19, 20) # nationalism

## best models
post.cov <- which(results.final$post.check == 1)
length(post.cov) ## 92 have postive definite covariance matrix
length(post.cov)/dim(results.final)[1] ## 10% of the models
good.results <- results.final[post.cov,]

head(good.results)
table(good.results$ndim)
## 1  2  3 
## 1 39 52 

##############
## Table 1
##############

## The best model: 3 dimensional
library(plyr)
best.results <- ddply(good.results, .(ndim), function(x){
  data.frame(x[which.min(x$chisq),])
})
(best.id<-which.min(good.results$chisq))
which.max(good.results$cfi)
which.max(good.results$tli)
which.min(good.results$rmsea)
print(good.results[best.id, ])
print(best.results[c(3,2,1),])
## p-values
1 - pchisq(301, 2)
1 - pchisq(2180, 3)

##############
## Figure 5
##############

pdf("graphs/cfa_chisq_dropQ6.pdf")
par(mar = c(4, 4, 1, 1))
plot(x = jitter(results.final$ndim[-post.cov]),
     y = results.final$chisq[-post.cov],
     xlim = c(0.5,7.5), ylim = c(64500,68100),
     col = "gray80", xlab = "", ylab = "", cex.axis = 1.2) 
points(jitter(results.final$ndim[setdiff(post.cov,best.results$id)]),
       results.final$chisq[setdiff(post.cov,best.results$id)],
       pch = 1, col = "gray30")
points(best.results$ndim, best.results$chisq,
       pch = 16, cex = 1.5, col = 1)
text(1, 67800, "Model C", cex = 1.2)
text(2, 65900, "Model B", cex = 1.2)
text(3, 65580, "Model A", cex = 1.2)
legend("topright", c("Valid","Invalid (neg. def. cov.)"),
       col = c("gray10","gray70"), cex = 1.3,
       pch = 1, bty = "n")
mtext("# Latent Factors", 1, 2.5, cex = 1.6)
mtext("Chi-squared", 2, 2.5, cex = 1.6)
graphics.off()

pdf("graphs/cfa_rmesa_dopQ6.pdf")
par(mar = c(4, 4, 1, 1))
plot(x = jitter(results.final$ndim[-post.cov]),
     y = results.final$rmsea[-post.cov],
     xlim = c(0.5,7.5), ylim = c(0.07385, 0.07560),
     col = "gray80", xlab = "", ylab = "", cex.axis = 1.2) 
points(jitter(results.final$ndim[setdiff(post.cov,best.results$id)]),
       results.final$rmsea[setdiff(post.cov,best.results$id)],
       pch = 1, col = "gray30")
points(best.results$ndim, best.results$rmsea,
       pch = 16, cex = 1.5, col = 1)
text(1, 0.0753, "Model C", cex = 1.2)
text(2, 0.07425, "Model B", cex = 1.2)
text(3, 0.07415, "Model A", cex = 1.2)
legend("topright", c("Valid","Invalid (neg. def. cov.)"),
       col = c("gray10","gray70"), cex = 1.3,
       pch = 1, bty = "n")
mtext("# Latent Factors", 1, 2.5, cex = 1.6)
mtext("RMSEA", 2, 2.5, cex = 1.6)
graphics.off()




